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ABSTRACT 

We investigate the formulation and application of element-level, element-independent 
error indicators. Our research culminates in the development of an error indicator formu- 
lation which is derived based on the projection of element deformation onto the intrinsic 
element displacement modes. The qualifier ‘element-level’ means that no information from 
adjacent elements is used for error estimation. This property is ideally suited for obtaining 
error values and driving adaptive mesh refinements on parallel computers where access to 
neighboring elements residing on different processors may incur significant overhead. In 
addition such estimators axe insensitive to the presence of physical interfaces and junctures. 
An error indicator qualifies as ‘element-independent’ when only visible quantities such as 
element stiffness and nodal displacements are used to quantify error. Error evaluation 
at the element level and element independence for the error indicator are highly desired 
properties for computing error in production-level finite element codes. Four element-level 
error indicators have been constructed. Two of the indicators are based on variational 
formulation of the element stiffness and axe element-dependent. Their derivations axe re- 
tained for developmental purposes. The second two indicators mimick and exceed the 
first two in performance but require no special formulation of the element stiffness and 
axe therefore element-independent. We describe an algorithm for element-splitting mesh 
refinement which we demonstrate for two dimensional plane stress problems. We discuss 
the parallelizing of substructures and adaptive mesh refinement. We demonstrate the final 
error indicator using two-dimensional plane-stress and three-dimensional shell problems. 



1. Introduction 


Automated solution techniques for the finite element method axe being implemented more 
and more frequently in the engineering sciences. As structures and structural models be- 
come increasingly complex, the need arises for techniques which require less interaction 
with the user while still returning accurate results. One particular advance in computa- 
tional techniques is the use of parallel processing for both solution procedures and adaptive 
mesh refinement. High speed, multi-processor computers can return solutions in a fraction 
of the time needed by serial computers. 

An important issue in adaptive mesh refinement is to establish the quality of the so- 
lution obtained in terms of discretization error. Many users of finite element analysis 
rely on convergence arguments to determine the accuracy or acceptability of a solution. 
In principle, meshes are refined until a convergent solution is obtained and that solution 
is generally accepted as correct. Another common procedure has been to quantify stress 
jumps between adjacent elements as a way to ascertain the quality of a displacement-based 
solution. These two approaches have obvious physical appeal to engineers. 

Quantifying discretization error is usually performed using a posteriori error indica- 
tors. An indicator is said to be element-independent if the method for computing an 
error value relies only on visible or derived quantities readily obtained from the element. 
For displacement-based finite elements, these quantities consist of stiffness, displacement, 
stresses, or strains. Element-independence is paramount in establishing error estimations 
for existing finite element codes in which element formulations may be considered black box 
computations. The present research has evolved from element-dependent error estimation 
to error analysis more suited to production-level codes. 

Error indicators are said to be element-level if they can be calculated without informa- 
tion from adjacent elements. This property is highly desirable in massively parallel compu- 
tations in which individual elements or element groups are stored on different processing 
units. Element-level estimation reduces processor co mm unication overhead, especially if 
neighboring elements are not stored on neighboring processors. 

Another important issue in adaptive mesh refinement is the choice of methods by which 
to refine the mesh (or parts of the mesh). Principally, there are four techniques by which 
to refine meshes: node relocation - r-adaptation; element splitting - h-adaptation; element 
rezoning - z-adaptation; and modification of the element polynomial order - p- adaptation. 
Consideration of the p- adaptation technique will be ommitted from this research as it 
requires tremendous changes in data structures and is therefore unsuitable for efficient 
parallel analysis. Two of the remaining three techniques have been implemented in previous 
research and we offer justification for the implementation of the third in the sequel. 

Decomposition of the finite element model into substructures is also a parallelizable 
process. Element groups are considered independently on individual processors where a 
solution step is performed. The solutions from several processors are then assembled and 
solved (usually via iterative methods) to create a complete solution for the structural 


2 


model. Hence, both the solution and mesh refinement phases of a complete analysis can 
be done in parallel. 

We defined three goals by which to guide the research: 

Goal 1. Development of an element-level error indicator capable of locating regions of 
error for single- and multiple-domain finite element meshes. The error indicator 
should not be affected by “jumps” due to material or geometric interfaces. 

Goal 2. Development of an element-independent error indicator which meets the objec- 
tives of Goal 1 and is capable of returning reasonable estimates for mesh dis- 
cretization energy error. 

Goal 3. Application of the indicator to solving realistic shell structure problems using 
parallel analysis. 

We present the progression from element-dependent error indicators to element-inde- 
pendent indicators. In the progression, the four error indicators are denoted sequentially 
as , C2, e 3 and € 4 . We demonstrate the utility of the two element-independent error 
indicators with regard to adaptive mesh refinement and assess their potential for measuring 
strain energy error. 

2. Historical Background 

The importance of developing reliable and economical finite element error estimates has 
been recognized by a large number of researchers since the late 1970s. Substantial progress 
in this direction has been made over the past decade. The following survey highlights key 
sources and advances in what is presently an active area of research. 

2.1. Bounding the Energy Error. 

The first attempt to establish a posteriori error bounds in problems of structural mechanics 
was made in the early 1960s by Fraeijs DeVeubeke 1 ’ 2 ; follow up work may be found in the 
Memorial Volume. 3 DeVeubeke advocated solving the structural problem twice, once with 
conforming elements and once with equilibrium elements, to obtain lower and upper bounds 
II_ and n + , respectively, to the actual potential energy II. The difference 


An = n + -n_ (2.1.1) 

may serve as an overall indicator of discretization error because An — * 0 as both meshes are 
refined. Bounds on pointwise errors can be derived through application of dummy forces 
and displacements and Castigliano’s theorem. Each such estimation requires, however, 
the complete solution of two linear problems. Despite its elegance, this approach did not 
attract attention from finite element developers and users. Its key drawbacks are: 

1. Two meshes have to be constructed. Because conforming and equilibrium elements 
have very different freedom configuration the two meshes do not have the same con- 
nectivity. 


3 



2. Equilibrium elements have never become popular because they are difficult to con- 
struct as well as expensive in terms of degrees of freedom and connectivity when used 
within a stiffness-based program. Similarly, purely conforming, exactly integrated 
elements are rarely used because they are too stiff. 

3. The computation of local error estimators, which are needed to drive mesh refinement, 
is cumbersome and expensive. Even if they had been economical to obtain, adapta- 
tion would involve simultaneously adapting two separate finite element models while 
satisfying locality constraints. This leads to significant implementation difficulties. 

2.2. Potential Energy Estimators. 

A second idea emerged in the early 1970: node relocation based on energy minimization. 

Consider a mesh of conforming finite elements with visible degrees of freedom collected in 

vector v. The potential energy 


n(v) = U(v) - W(v) = ±v r Kv - p T v, (2.2.1) 

is considered also a function of the node locations: II(v,x). In the linear case, II = 
^v r K(x)v — p r (x)v. For fixed x, the solution v = K -1 p gives II(x) = p(x)K(x) _1 p(x). 
This can be minimized with respect to x to get an “optimal mesh” with the given topology 
and freedom configuration. This approach is now called r-adaptive mesh refinement, or 
the r-adaptation for short. 

2.3. Local Estimators. 

The idea leading to the selfadaptive h-method, or h-adaptation for short, was appar- 
ently first proposed in two papers by Babuska 4 and Sewell 5 that appeared on the same 
MAFELAP Proceedings volume. The article by Babuska was, however, more specifically 
oriented to applications and thus exerted bigger influence. These first error estimators for 
h-adaptation were of residual type. For the prediction phase many error indicators and 
error estimators have been developed over the past 15 years. The most successful ones are 
based on measuring discontinuities of derived field quantities such as stresses or strains 
across interelement boundaries. Other estimators Me based on interpolating residual es- 
timates considering element patches. In 1986 Zienkiewicz and Zhu 6 proposed an error 
estimator based on the smoothing of the stress jumps at element nodes. More recently, in 
1992, Zienkiewicz and Zhu 7 proposed an improved error estimator with superconvergent 
characteristics. This has been widely adopted. 

2.3.1. Zienkiewicz- Zhu Superconvergent Patch Recovery 

We now briefly review the Zienkiewicz- Zhu (Z 2 ) superconvergent patch recovery for error 
estimation. 7 The basic measure used in Z 2 error estimation is an energy norm of the 
difference between a smoothed solution over an element patch and the derived values from 
the element. A patch is defined as a group of elements surrounding a node. In essence, 
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smoothing a solution over a patch effectively elevates the polynomial order by one. We 
confine our discussion to simple problems of linear hyperelasticity. 

Define the smooth solution a* as 


( 7 * = Pa. (2.3.1) 

where P is a polynomial basis representing the assumed displacement field of the element 
and a is a vector of unknown parameters. Minimizing the function 


F( a) = ^2 { a ( x iiyi) ~ P(*;>2/;) a ) 2 (2.3.2) 

1=1 


implies that 


a = A 1 b 


(2.3.3) 


where 

n 

A = ^P T (xi,T/i)P(x, -,</,), 

t=l 

n 

b = ^2P T (xi,yi)<r(xi,yi). 

1=1 


(2.3.4) 


The number of sampling points, n = mk, is equal to the number m of elements surrounding 
a node times the number of sampling points k per element. 

For two and higher dimenisonal problems, the Z 2 error estimate e z on the element 
level is given as 

||e z2 ||, = / (<T-<r’) T (<r - <r‘)i{l t (2.3.5) 

where <r and <r* are vectors of element- derived and smoothed stresses, respectively and f2 e 
is the volume given of the element. Averaged nodal values used for error estimation are 
given by 

n«*’iu = {eoi^’ii.)’}’- (2 - 3 - 6) 

k «=i 


By now there is a large and rapidly growing literature on FEM local error estimators. 
The surveys by Babuska 8 , Zienkiewicz e< a/. 9 ’ 10,11,12 and the finite element book by Szabo 
and Babuska 13 may be recommended. In summary, it can be stated that present estima- 
tors work well for conforming elements, homogeneous domains and continuum-type linear 
problems. The reliance on jump conditions or the smoothing of the jumps brings difficul- 
ties, however, when such jumps are a natural part of the problem as in junctures, material 
interfaces, localization bands or shocks. Moreover, conforming elements are not necessar- 
ily the best performers in structural problems. We now investigate the development and 
application of four element-level error indicators. 
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3. Parameterized-Functional Error Indicators 

Our original research of element-level error estimation is based on extraction of higher 
order element energy from high-performance elements. As described below, these ele- 
ments are based on parameterized functionals and the development of two error indicators 
through specialized formulation has lead to Parameterized-Functional (PF) Error Esti- 
mation. While these particular error indicators require element-dependent formulation, 
their derivations have significant implications for the present research and are retained for 
developmental purposes. 


3.1. Background 


Felippa and Militello 14 ’ 15 introduced a general parametrized form of the stress- strain- 
displacement functional of linear hyperelasticity. A general expression for the functional 
is: 


II = U a fo - P = U( u,a,e,a,P,i) - P, (3.1.1) 


where U is the generalized strain energy stored in the body volume, P is a forcing potential 
that embodies all other actions such as loading and boundary conditions, a, 0 and 7 are 
free parameters further discussed below, and u, cr and e are the varied displacements, 
stress, and strain fields, respectively. (A superposed tilde expresses that the symbol is 
subject to independent variation.) The variation 511 = 0 generates all governing equations 
of elasticity for arbitrary a, 0 and 7 . All elasticity functionals with varied displacements 
can be obtained by appropriate specialization of these parameters. 

As to the forcing potential P, three forms labelled as P c , P d and P* have been studied 
in conjunction with parametrized functionals . 16,17 ’ 18 ’ 19 * 20 P c is the conventional forcing 
potential, which is used in non-hybrid finite element discretizations. Potentials P d and 
P‘, called displacement-generalized and traction-generalized, respectively, are useful in 
the construction of high-performance hybrid elements with additional boundary fields. Of 
these the displacement-generalized potential P d has proven to be the most useful one in 
the construction of high-performance elements . 21 This potential depends on three varied 
fields 


pd 


P rf (u,d-,d), 


(3.1.2) 


where d is a boundary displacement field defined only at element interfaces. A common 
property of P c , P d and P‘ is that they do not directly depend on the free parameters. 

Now suppose that (3.1.1) is used to construct a finite element discretization, while 
keeping one or more parameters free. (Practically this means that such parameters are 
kept as arguments of the element stiffness subroutines.) The discrete approximants u, cr 
and e obtained for a given mesh will then be generally functions of the free parameter(s). 
An important property of (3.1.1) is that U takes the same value if the varied fields u, 
ff, e are set to the exact solution fields u, <r, e, independently of a, 0 and 7 . Hence, we 
may expect that the difference between two values of U obtained for two different sets 
of parameters may provide a indicator on how far we are from the converged solution. 
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That is, the difference may be adopted as an error indicator. The most interesting feature 
of this indicator is that it naturally provides an element-level indicator, as the following 
development shows. 

To tackle the general case first, suppose that the two sets of parameters are (oti, /?i, 
7 j) and (a 0 , /?oi 7 o)- The corresponding approximate values for displacements, strains 
and stresses obtained with a given mesh are (ui, 07 , ei) and (uo, 0 o, eo), respectively. 
Denote by = U e ( u x , 07 , § 1 , 01 ,^ 1 , 71 ) and U% = f7 e (u o ,0- o ,eo,ao,^o,7o) the value of 
the generalized strain energy evaluated over the element of that mesh. Then the PF 
element error indicator is defined as the difference 


A = \UZ~UZl (3.1.3) 

where p denotes a parameterized-functional error indicator. This definition apparently 
requires that the problem be solved twice for a given mesh. It will be seen, however, that 
consideration of the structure of the stiffness matrices and their dependence on the free 
parameters allows the use of only one solution. 

Let K e be the stiffness matrix of an ANDES 22 element. This is a variant of the 
Assumed Natural Strain (ANS) formulation, a name coined by Park and Stanley . 23 Let v e 
be the visible element degrees of freedom (those degrees of freedom in common with other 
elements, also called the connectors) and p the corresponding element node forces. Then 
the element stiffness equations decompose as 

K e v e = (KJ + Kfc) v e = p e . (3.1.4) 

K‘ b and K\ are called the basic and higher order stiffness matrices, respectively. The basic 
stiffness matrix, which is usually rank deficient, is constructed for consistency. The higher 
order stiffness matrix is constructed for stability and (in more recent work) accuracy. A 
decomposition of this nature, which also holds at the assembly level, was first obtained by 
Bergan and Nygard 24 in their derivation of the Free Formulation. 

3.2. The First PF Error Indicator 

With a view towards reducing the need for two solutions, the first PF indicator, PF 1 , is 
chosen as the difference between the element internal energies for aj = a and c*o : 



(3.2.1) 


Next we assume, without proof, that for two different parametrizations and for a ‘good’ 
mesh the following property holds: 0 * a ~ cro, u a ~ uo, and e Q ~ eo; but that e a ^ ^ 

e “> e 0 / ej / e#. In other words, corresponding fields in two different parametrizations 
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converge toward each other faster than do different fields for each parametrization. If this 
assumption holds we can approximate the error indicator (3.2.1) as 

e\ = \a f {&- o*)Z(e 9 - e) tt dV. (3.2.2) 

Jv 

Through algebraic manipulations this expression can be transformed 21 to 

= i(v‘ ) r KJv*, (3.2.3) 

where v* denote the element degrees of freedom calculated for the a solution. Physically 
(3.2.3) is the higher order energy (HOE) absorbed by the element. Since for ANDES 
elements (as well as for FF elements) the higher order stiffness matrix K£ is always available 
in separate form, this indicator is readily calculated in an element by element manner. 

It should be noticed that the PF1 error indicator is related to that heuristically proposed 
by Melosh and Marcal 27 in the context of the SED (Strain Energy Density) method. 

3.2.1. Important Properties of the PF1 Error Indicator. 

The error indicator inherits the properties of the higher order stiffness. By construction 
the higher order stiffness verifies 

KM = o, (3.2.4) 

where vj are the nodal displacements associated with a rigid body motion or a constant 
strain state. Actually Vj expands the null space of K*. Any other v e can be decomposed 
a sv* = Vj + vJ, where v£ is the projection of v e on the null space of K e h . 

From these considerations we conclude that (3.2.3) automatically “filters out” the con- 
tributions to the displacement field from rigid body motions and constant strain states. 
Thus, in problems whose analytical solution consists of uniform strain states the error 
indicator vanishes over each element, in accordance with the fact that any mesh should 
solve those cases exactly. This statement may be extended to cases solved by constant 
strain states that jump at material or geometric interfaces, if such interfaces axe modeled 
exactly by the finite element mesh. Other error indicators in present use do not verify this 
property. 

3.2.2. Performance of PF1 

The PF1 indicator has been tested extensively in numerical and adaptive mesh refinement 
studies 28 ’ 29 and is capable of capturing regions containing high strain gradients. However, 
the error measured by PF1 does not converge to the true strain energy error and is element- 
formulation dependent. Thus the requirements of Goal 2 are not met. With a view towards 
constructing an estimator which will return the true energy error, we turn to the next PF 
error indicator, PF2. 
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3.3. The Second PF Error Indicator 


As previously stated, the PF1 indicator has the ability to locate high strain gradients 
near junctures and singular points. The PF1 indicator is useful for error indication and 
adaptive mesh refinement in problems which exhibit these characteristics but fails to return 
an accurate measure of strain energy error. We now describe the second PF indicator, PF2. 
PF2 partially fufills the requirements of Goal 2 in that a true measure of the element strain 
energy error is returned for some displacement fields. While this indicator has remains 
undeveloped (primarily due to its inherent formulation dependency), its derivation and 
application in numerical experiments have lead to the most recent form of the element- 
independent error indicators. 

3.3.1. Determining Error by Energy Balance. 

In order to estimate the true energy error we need to know precisely how an element 
behaves in a particular strain field. We now examine behavior of the three node Extended 
Free Formulation (EFF) high-performance triangle element 30 in a pure bending field. To 
determine the energy error for an EFF element we begin by applying a known displacement 
mode (in this case, a pure bending field) to an element patch as shown in Fig. 1. The 
analytical strain energy for this element patch is 

UdV (3.3.1) 

where U is the strain energy density for the given displacement field. Strain energy for the 
finite element patch is given as usual by 

U fe = |v eT K e v*. (3.3.2) 



K xx = 10 


Element Aspect: H 

Fig. 1. Element patch of two EFF membrane triangle elements: 
Pure bending field. 
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The true strain energy error is defined as the difference between the analytical and 
finite element patch strain energies 


ttr = \U* n -U U \. (3.3.3) 

The element stiffness matrix K* for the EFF triangle is decomposed as 

K'=Kt(a b ) + p b K e h (3.3.4) 

where a b is a variable scaling factor for the force-lumping scheme used in computing Kj 
and fl b is an overall scaling factor for K*. We wish to obtain a new stiffness matrix K ,e 
such that the element error is also given by 

e 2 = §( v fc) TK 'fc v fc> (3-3.5) 


Symbolic solution for this element patch (setting ej r = e|) in pure bending yields two 
unknowns in two equations as coefficients for the element height H and height cubed H 3 . 
From these equations we return a relationship from which a' and /?' may be computed 
given a b and 

2 a b — 3 


a = 


/?' = 


5 - 4 a b + 20 b + 4t/ 2 


(3.3.6) 


The element strain energy error is then determined using (3.3.3). 
3.3.2. Behavior of the PF2 Error Indicator 


For the two-element patch in pure bending we are able to return a true measure of the 
strain energy error. However, numerical studies reveal that PF2 is highly sensitive to 
components of displacement modes other than pure bending. For example, Fig. 2 shows 
the effect of rotation of the third node beyond the pure bending field rotation. The rotation 
of node three is swept between ±10% of the pure bending rotation. Note that for Poisson’s 
ratio v = 0, aj = 3/2 and /3 b = 1, we obtain a' = 0 (constant strain triangle stiffness) and 

P = 1 / 2 - 
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Fig. 2. Effect of Node 2 rotation beyond pure bending rotation: 
a' and ft values. 

Clearly, the additional drilling rotation has a profound effect on the a' and ft values. 

Similar effects occur for other displacement modes applied to the element patch. 

S.S.S. Toward an Element- Independent Error Indicator 

From the PF1 and PF2 error indicators, we have established several important points 

regarding element behavior and error estimation: 

1. All strain energy error is contained in the higher order element displacement modes. 
Therefore, filtering out basic (i.e., constant strain and rigid body) modes is desirable 
for computing an error measure. No contribution from basic modes also allows for error 
estimation despite “jump” conditions due to material and geometric discontinuities. 

2. Error indicators are apparently sensitive to various displacement mode configurations 
of the element. This point requires us to quantify contributions to strain energy error 
from (higher order) modes experienced by the element. 

3. High strain gradients are captured by both PF error indicators enabling these indi- 
cators to be used for structural problems in which analysis near singular points is 
desirable (e.g., linear-elastic fracture problems). 
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4. Element-Independent Error Indicators 

The PF indicators provide an excellent base upon which to build error estimation ideas. 
The primary shortcomings of PF indicators and estimators are, of course, their element- 
formulation dependencies. In order to extend our error estimation ideas to standard finite 
elements, we introduce Element Independent (El) error estimation. We have again con- 
structed two error indicators in an attempt to mimick and then improve upon the PF 
indicators. 


4.1. The El Error Indicator 

First, we wish to construct an El error indicator that is capable of matching the per- 
formance of PF1 (i.e., one that will capture high strain gradients near junctures and 
singularities). Realization of this objective satisfies the element-level requirement stated 
in Goal 1. Development of this indicator, Ell, provides a step for our proceeding to meet 
the requirements in Goal 2. 

4.1.1. Development of Ell 

As previously shown, a property of the PF indicator is that strain energy due to lower strain 
modes is filtered out. This means that higher order strain energies are used to indicate 
element error. The question is then raised, “Can higher order strain energy from any finite 
element be used as an error indicator?” This is possible since most finite elements will 
produce higher strain modes, a method should be developed by which higher order energies 
for those elements could be extracted. With this thought in mind, we construct the Ell 
indicator with characteristics similar to those of PF1. 

Let us examine a decomposition of the total strain energy 

U = |(vfKv‘. (4.1.1) 

Splitting the displacement vector v into basic and higher order modes we can re-write the 
total strain energy as 


i(v e ) r Kv e = i(vg) r Kvg + (vg) T Kv£ + ±(v£) t Kv£. (4.1.2) 

The basic strain energy is 

U„ = i(yi) T KvJ. (4.1.3) 

Note that (4.1.2) is the energy equivalent of the total strain energy min us the HOE from 
the PF indicator. This is because vj expands the null space of the higher order basis of 
K. We define the augmented element HOE as 


Vi = iK) r KvJ 

= i(v'f Kv' - (1 (vJ) t KvJ +(vJ) i ’Kv;) 


(4.1.4) 
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where the “+” indicates an additional energy term. The algebraic difference between 
and the previous HOE Uh (used as the PF1 indicator e\) is 

of - U h = {(vJ) T (K;)-vt (4.1.5) 

where (K£)* would be the basic element stiffness if such could be found. 

Numerical experimentation shows that the term 7 (v^) t (Kj)*v^ makes U£ too con- 
servative to be used for reasonable error indication. However, we can use as part of a 
scaling factor for the total strain energy to return reasonable values for error indication in 
regions of high strain gradients and in the vicinities of singularities. 

We now define our Ell indicator as 



(<f Kv£ 

(vS) t Kv£ + K) T Kv‘ 


(v*) r Kv‘ 


— toti 


(4.1.6) 


where a& is a scaling factor for the total strain energy Utot and the superscript i indicates 
element-independence. Note that the element stiffness K remains intact for the error 
computation and requires no decomposition. In addition, the numerator of in (4.1.6) 
contains only energies from higher order modes, meeting the requirement that “jump” 
conditions in constant strains will not affect the error measure. 

The principal advantage of Ell is that it can be used for error indication in standard 
finite elements, provided v* and v& are readily found. We now turn our attention toward 
decomposition of the displacement vector and computation of the necessary energy terms 
to construct Ell. 

4.1.2. Obtaining Higher- Order Energy using Projectors 

For simple planar elements (e.g., 3-node triangle), obtaining the higher order displace- 
ments v£ is relatively straightforward. To do this, we use a projector which filters out 
displacements for rigid body and constant strain modes. The projector, developed by 
Bj0m Haugen 31 , is now defined. Recall the decomposition, 


V e =v£+v£. (4.1.7) 

We write vj as a linear combination of the rigid body and constant strain (basic) mode 
vectors 

vl = Ra (4.1.8) 

where R is a basis for v£ and a is a vector of coefficients such that 


v e = Ra+ vj. 


(4.1.9) 


We require 
such that 


that the higher order displacement vector v£ 


R r vt = 0. 


be orthogonal to the basic modes 


(4.1.10) 
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Premultiplying v* by the transpose of our basis matrix R and solving for a yields 


a = (R r R)" 1 Rv*. 

Upon rearrangement using (4.1.11) and (4.1.9) , we obtain 

v% = (I - R(R T R) _ 1 R T )v e 
= Pv e . 


(4.1.11) 


(4.1.12) 


Note that P 2 = P, which characterizes a projector matrix. The basic energy and aug- 
mented HOE are then computed using (4.1.3) and (4.1.4), respectively. 

The projector algorithms for 3-node triangle elements have been implemented and 
tested to ensure that they extract the higher order displacement modes. They are now 
in use for computation of the Ell indicator. Projectors for 3-node, 4-node and most 
planar elements are readily constructed with a second projector sometimes necessary, as 
in the case for warped 4-node elements. The projector method can also be extended to 
three-dimensional elements such as 8-node bricks. For more sophisticated elements with 
non-colinear nodes (e.g., doubly-curved shell elements), we will need other methods of 
determining error energy. 

4-1.3. Obtaining Higher- Order Error Energy using Strain Averaging. 

For higher order elements, the computation of error energy becomes somewhat more com- 
plex. In filtering out the basic modes by use of projectors, we are in effect subtracting 
an average strain energy from the total energy. To extend this concept to higher order 
elements, we write the basic energy in its integral form 

Ub = k [ & T edV. 

Jv 

and also write the augmented HOE in its integral form 

U h = \ [ (<r - d-) T (e - e)dV 

Jv 

where e and <r are the total strain and stress and e and & are the averaged strain and 
stress, respectively. Define the averaging strain-displacement matrix as 

B = i/ v B(e,,)dV 

with 

v= f [J«,r?)| dtir, (4.1.16) 

Jv 


(4.1.15) 


(4.1.14) 


(4.1.13) 
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and £ and rj axe the element natural coordinates and J(£, tj) is the transformation Jacobian. 
The higher order strain and stress vectors are then 


e h = (e — e) = (B — B)v, 
v h = (cr - &) = Ce h 


(4.1.17) 


with C being the element constitutive matrix. Our El error indicator is then defined as 

« - (41 - 18) 

As presented, this method requires knowledge of the element strain-displacement matrix 
B. This detracts somewhat from the idea of a truly element-independent error indicator. 
However, in lieu of using the actual B matrix, we may assume a strain-displacement 
relationship B* provided we have a priori knowlege of the general element displacement 
capabilities. Our error indication may then be obtained by substituting B* into (4.1.15) 
and (4.1.16) and computing t\ using (4.1.18). 

4.2. The Deformation Mode Error Indicator 

We have seen from the development of the PF1, PF2 and Ell indicators that strain energy 
error occurs in higher deformation modes and that the error indicators appear to be sen- 
sitive to how much of the higher modes experienced by the finite element. Therefore, we 
propose that quantifying the strain energy error necessitates determining amounts of error 
contained in the higher modes and the extent to which each of the modes is represented 
for a fini te element in a given displacement field. We introduce a new error indicator based 
on projection of the element deformation onto the intrinsic element displacement modes 
which we call the Deformation Mode (DM) Error Indicator. 

The DM error indicator offers promising results toward realizing both Goals 1 and 2. At 
present, DM has the capability of returning true energy errors in some displacement fields 
and a reasonable representation of error in others. We begin with the derivation of the 
indicator for planar two-dimenisonal elements and evaluate performance of the indicator 
with example problems in the sequel. 

4-2.1. Derivation of the Deformation Mode Error Estimator. 

A basis for the displacement modes of a finite element is given by the eigenvectors of the 
element stiffness matrix 

¥ = eig(K). (4.2.1) 

An alternate basis is determined considering the known-mode displacements, which have 
a physical appeal with regard to a priori knowledge about continuum mechanics for those 
displacements. We define this second set of eigenvectors as 

* = eig(K) (4.2.2) 
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where A such that 


diag(A) = * T K# 


(4.2.3) 


are the eigenvalues of K as usual. 

For the present, we assume that the displacement field within an element domain is 
comprised of a Unear combination of the known displacement modes of the finite element. 
Obviously, the finite element approximation to the displacement field is comprised of these 
modes. 

The analytical strain energy for any of these modes is readily computed by integrating 
the strain energy density Uk for that mode over the element domain 

UL = f UkJV, (4.2.4) 

Jv 


where k = 1, . . . , n specifies which of the n modes is being considered. The strain energy 
for the finite element is 


Ufe = f(v*) r Kv* 

and the strain energy error for the k th mode is given as 

(4.2.5) 

= U k an - U%. 

(4.2.6) 


Choosing an acceptable displacement field for the analytical field and requiring that 
the finite element match that field at the nodal points, we compute the error energy for the 
element in the prescribed mode using (4.2.4). We now have a scalable error value for the 
finite element in each of its modes. These error values are then coUected into the modal 

error vector 

m 1 

m 2 
m* 

We propose that the apparent energy error for the finite element is a superposition of the 
amounts of error contained within each of the displacement modes. 

To determine the amount of error attributable to each mode, the finite element displace- 
ment vector is decomposed into its modal components using the inverse of the eigenvectors 
of the stiffness matrix 

<*m = $ -1 v. (4.2.8) 

The modal component vector a m represents the magnitude of each of the modes seen in 
the element displacement field. The finite element strain energy is then written 

u ft = K&at m . (4.2.9) 
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Recalling the assumption that the analytical displacement fields can be represented by 
the eigenvectors $ of K, and that the element nodal displacements can be decomposed 
into their respective modes, we have a means by which to compute an error measure for 
the element. Because the eigenvectors $ diagonalize K, we can write the error for the k th 
mode as 

e* = crfm*. (4.2.10) 

The total DM error e\ is then obtained by 

= a\m x + a\m 2 + . . . + a 2 n m n . (4.2.11) 

where a, and m,- are the modal component and modal error vectors, respectively. 

4-2.2. Application of DM to Shell Elements. 

Obviously, some finite elements can represent displacement fields for which analytical com- 
plements axe not easily found. This problem is resolved by “borrowing” surface modeling 
techniques from computer- aided design and manufacture (CAD and CAM). 

Here, we define a 4-node element in parametric form using a tensor-product Ferguson 
surface patch (FSP). 32 The FSP is used els it relies on corner conditions which are derived 
readily from nodal values of the finite element. Fig. 3 shows a FSP with the corner 
conditions depicted in vector form (subscripts indicate derivative quatities and superscripts 
denote corner points). 



Fig. 3. Ferguson Surface Patch. 
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We define the FSP equation as 


r,(a,t) = SCQC t T t 


(4.2.12) 


where 


S= [1 S S 2 S 3 ] 0 < 3 < 1, 

T = [ 1 t t 2 t 3 ] 0 < t < 1, 

r i o o oi 


-3 3 -2 -1 ’ 
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dadt l°° 

dadt I 01 

dr | 

ar | 

a 3 r | 

a 3 r 1 

3711° 

aTlii 

3I3tlio 

dsdt 


(4.2.13) 


Here, P and r-derivatives denote the ( x , y, z ) coordinate triples and comer configurations 
for the FSP, respectively. Subscripts 00, 01, 10 and 11 specify the comer points for which 
each value is given. 

The eigenvector modes of standard shell finite elements contain no information about 
the cross derivatives ( “twists”) at the node points which are necessary to complete the lower 
right 2x2 portion of Q in (4.2.13). To obtain a reasonable model for the finite element 
displacement modes, we need to estimate these twists somehow. As we wish to compute 
energy error for each of the finite element modes, we will use an energy- minimizing twist 
estimator based on a variational approach. A survey of other twist-estimation procedures 
is given by Farm. 33 


The functional 



(4.2.14) 


is given by Nowacki, Reese 34 and Walter 35 as a fairness criterion for surfaces in engineering. 
The minimum and maximum normal section curvatures are /e min and « mar , respectively. 
This functioned is used because it is a measure of the strain energy of flexure and torsion 
for a thin rectanguleir plate under small deflections. Hagen and Schulze have used this 
functioned for twist estimation of patches with orthogonal boundary curves. 36 Farin has 
extended it to include nonorthogoned patch boundaries giving an estimate for the normal 
component of twist at a point as 37 


where 


T d 2 r _ h 
n dsdt g 


(4.2.15) 


dr T dr ,dr T dr T d 2 r dr T dr T d 2 r, 

~drdt^~dTd^ n d^ + ~dTdt n ds 2 ' 


(4.2.16) 
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(4.2.17) 


_ dr T dr dr T dr , dr T dr . 2 

^ ds ds dt dt ds dt 
The tangential components of the twist vector are presumed small and set to zero. 

From the preceeding formulation we construct element “surfaces” for the deformed and 
undeformed shapes of the element, r(s,t) and r°(s, t), respectively. The deformation field 
u(s,t) is defined as 

u(s,t) = r(s, t) — r°(s,t) (4.2.18) 

from which strains and strain energies can be derived. The analytical and finite element 
strain energies for a given mode are computed using (numerical integration of) (38) and 
(39), respectively. 

4-2.3. Additional Energy Terms. 

The true displacement field for a mesh region may of course contain displacements not 
represented by the displacement modes of the finite element. Including the strain energy 
from these displacements we rewrite (45) as 

€t = OjCi + 02^2 + • • • + + Ua n, (4.2.19) 

where U an is the strain energy not included in the finite element approximation. As- 
sumedly, the additional strain energy U an will come from higher order or cross- term dis- 
placement modes and may or may not decrease with mesh refinement. 
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5. Example Problems 

We now demonstrate the applicability of the two element-independent error indicators. For 
the El indicator we examine a simple linear elastic fracture problem and demonstrate the 
h-adaptation refinement. Additionally, we show the performance of the DM indicator in 
membrane- and bending-dominated problems. We will see that the El indicator performs 
well for finding regions of high strain gradients and that the DM indicator qualitatively 
matches the Z 2 error estimator. 

5.1. Example of Ell: An Application to h-adaptation Mesh Refinement. 

We have shown how higher order energy can be extracted from a finite element to produce 
an indication of high strain gradients. We demonstrate the following h-adaptation algo- 
rithm in the context of a simple linear fracture problem given in the Example Problems 
section. 

In previous work, 28 ’ 29 we describe and demonstrate the r-adaptive (node relocation) 
mesh refinement method. As an accompaniment to our research involving error indicators 
is used in element-splitting of three-node triangle elements. 

Briefly, the logic used for h-adaptation is as follows: 

1. Flag a percentage of the elements which exhibit the highest error for “splitting.” 

2. Inject nodes on the edges of the elements flagged for splitting. Assign surface param- 
eters and boundary conditions to the new nodes. 

3. Determine the number of new elements required and insert them accordingly. This is 
performed in passes to avoid creating “island” elements - unsplit elements surrounded 
by split elements - which create incompatibility. 

4. Recompute error and follow-up with additional refinements as necessary to achieve a 
desired tolerance or convergence. 

The h-adaptation method is essentially the division of single elements to multiple ele- 
ments. This effect may prove useful in augmenting error measures as the solution converges. 

Fig. 4 shows a standard linear-elastic fracture problem of a simply supported beam 
with a vertical crack. Figs. 5 and 6 show the original mesh and the meshes resulting from 
two h-adaptations. Mesh refinement occurs as expected near the cracktip, applied load 
and boundary conditions. 

The Ell error indicator performs well in determining the areas of high strain gradients 
for refinement. We conclude that while the Ell indicator will not return proper strain 
energy error distribution, it can be successfully used in problems involving singularities. 
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Fig. 4. Linear Elastic Fracture Problem. 

Original mesh has 160 3-node triangle elements. 



Fig. 5. Linear Elastic Fracture Problem. 
Mesh after one h-adaptation. 



Fig. 6. Linear Elastic Fracture Problem. 
Mesh after two /i-adaptations. 
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5.2. Examples of DM Using 4-node Isoparametric Elements 

We now demonstrate the DM indicator for three plane stress problems using four-node 
isoparametric plane stress quadrilateral elements. In lieu of deriving the higher order mem- 
brane displacement fields, via the FSP, we will assign pure bending fields (i.e., quadratic 
fields) based on the two higher mode eigenvectors for this element. 

5.2.1. DM Error Indicator: Membrane Cantilever with Applied End Rotation 

Our first problem used to assess the DM indicator is a membrane cantilever, built-in on 
one end with an applied rotation to the opposite end. Four meshes were used to model 
this problem. The applied rotation creates a constant curvature for the analytical solution. 
Fig. 7 shows one of the meshes (4 x 40 plane stress four-node quadrilateral elements) used 
to test DM in a pure bending field. Aspect of the cantilever is 10 : 1. 
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0 = 0.01 
E = 1.0 
v =0.0 
t= 1.0 


Fig. 7. Membrane Cantilever Problem (2 X 20) elements: 
Built-in end with applied end rotation. 

For all elements: Uf e = 0.011719, U an = 0.010417, 
e an = £/e = 0.001302. 


DM returns the true strain energy error for each of the meshes tested (grids of 1 x 10, 
2 x 20, 4 x 40, and 8 x 80 elements). These tests, and those conducted on single elements 
in the development of DM, verify the ability of DM to decompose the displacement field 
v into its componenent modes. 

5.2.2. DM Error Indicator: Membrane Cantilever with Applied End Shear 

Our second example for the DM indicator is a membrane cantilever, shown in Fig. 8. The 
analytical displacement solution is cubic leading to a quadratic strain field. We can assess 
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the effect of U of (40) on the error estimate as we uniformly refine the mesh as the strain 
energy due to shear is not modeled by the four node quad. 

We start with a coarsely-meshed model with 10 elements. Fig. 9 shows that the DM 
indicator slightly overestimates the true error. Next we refine the mesh by halving the 
element size (40 elements). Now DM begins to underestimate the true error as shown in 
Fig. 10. The finite element strain energy is within about 12% near the root end and the 
computed free end tip displacement is about 90% of the analytical value. 

Continuing to uniformly refine the mesh we see that the error indicator begins to greatly 
underestimate the analytical error. In Fig. 11 we see that the high energy errors near the 
outer edges of the cantilever at the root end axe not captured by the DM indicator. Because 
of the uni form mesh refinement, the stiffness K e for all elements is the same from mesh to 
mesh. Thus we conclude that indeed U of (40) exists and is comprised of strain energies 
not captured by the DM indicator. 

However, it is also worth noting in the regions of apparently high error, that the true 
error has decreased to less than 4% of the total strain energy and that the tip displacement 
is 97% of the analytical value. DM also consistently predicts higher error (though under- 
estimated) toward the root of the cantilever. This means that mesh refinement would still 
tend to refine the mesh in the regions of highest error. 




P=0.1 
E= 1.0 
v =0.0 

t= 1.0 


Fig. 8. Membrane Cantilever Problem: 

Built-in end with applied end shear. 
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Fig. 9. Membrane cantilever with end shear, (1 X 10) elements: 

Comparison of true strain energy error to DM indicator error. 


True Strain Energy Error 
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Fig. 10. Membrane cantilever with end shear, (2 X 20) elements: 

Comparison of true strain energy error to DM indicator error. 
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Fig. 11. Membrane cantilever with end shear, (4 X 40) elements: 

Comparison of true strain energy error to DM indicator error. 




5.2.S. DM Error Indicator: Sharp Strain Gradients 

We now observe the performance of the DM indicator for capturing high errors in regions 
of sharp strain gradients. This problem was proposed by Oden et al. 3S and used by 
Zienkiewicz and Zhu 7 for testing their error estimator. 

The proposed problem is given as 


tt = a;(l - x)y(l - y)tan *(<*(£ -£o)) 
with recommended values of 

_ Qc + y) 

V2 

£o = 0.8 
a = 20. 

The governing equation of the problem is 

— Au = / 


(5.2.1) 


(5.2.2) 


(5.2.3) 


with boundary conditions u = 0 on dQ where Q is a unit square domain on (0,1) x (0,1). 
Figs. 12(a) and 12(b) show the contours of du/dx and du/dy, indicative of the strain 
gradients in x and y. 
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Figs. 13, 14 and 15 show the performance of the DM indicator with respect to capturing 
sharp strain gradients while uniformly refining the mesh: 5 x 5, 20 x 20, and 50 x 50 isopara- 
metric quadrilateral elements, repectively. Referring to Fig. 12, we see that increased error 
is predicted in the regions accompanying the high gradient. Addtionally, visual compari- 
son with the exact solution and patch-recovery error estimator of Zienkiewicz and Zhu 7 is 
favorable. 



Fig. 13. Sharp strain gradients (Oden): 
DM error indicator 
(10 X 10) elements. 



Fig. 14. Sharp strain gradients (Oden): 
DM error indicator 
(20 X 20) elements. 
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Fig. 15. Sharp strain gradients (Oden): 

DM error indicator 
(50 X 50) elements. 

5.2.4 • Exam-pie of DM Using 4-node ANS Shell Elements 

We now demonstrate DM using the benchmark cylinder-plate juncture (CPJ) problem 
developed in previous work. 29 Fig. 16 gives geometric, material, load and B.C. data for 
the CPJ problem. We compare our DM indicator with a refined model tested at the 
Structures Laboratory at Lockheed PARL. 39 All analyses of the CPJ problem are modeled 
using quarter symmetry. 

6.667 w ‘ 



Fig. 16. Benchmark problem featuring cylinder-plate junctures. 

The cylindrical shell is supported by end diaphrams, with 
fixed derees of freedom as shown. 
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Plate 1 shows the relative error distribution determined by the Lockheed testbed code 
for a mesh of 320 16-node ANS shell elements. For all plates, the error ranges blue to 
red from low error to high error, respectively. The Lockheed solution is based upon the 
Zienkiewicz-Zhu 6 error estimator using special mesh partitioning to remove the effects 
of the juncture. Plates 2 and 3 show the performance of the DM indicator for meshes 
of 320 and 1056 4-node Assumed Natural Strain - Equilibrium Constrained (ANS-EC) 
elements. 40 The ANS-EC 4-node element belongs to the class of C° plate/shell elements. 
Higher-order contributions to the element stiffness axe given from a strain-displacement 
matrix which is not derived from the standard Lagrange shape functions, but is instead 
derived from the imposition of kinematic and equilibrium constraints. 

We see from Plates 2 and 3 that the DM indicator qualitatively matches the results 
from the Lockheed error estimation shown in Plate 1. The mesh in Plate 2 uses the same 
number of elements as used for the Lockheed model in Plate 1. Plate 3 shows a similar 
error distribution along the juncture. Note that in all cases, the error is primarily confined 
to the single rows of elements directly on either side of the CPJ. To date, the DM indicator 
for shell elements has not been tested quantitatively to establish its potential for true error 
estimation. 
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Plate 1. Benchmark problem featuring cylinder-plate junctures. 

Lockheed PARL relative error estimation for 320 16-node ANS elements 
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Plate 2. Benchmark problem featuring cylinder-plate junctures. 

DM error indicator for 320 4-node ANS-EC elements. 
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Plate 3. Benchmark problem featuring cylinder-plate junctures. 

DM error indicator for 1056 4-node ANS-EC elements. 
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6. Summary 

We have developed an element-level, element-independent error indicator known herein 
as the DM error indicator. This new error indicator is founded upon ideas stemming 
from two previous error indicators based on parameterized variational principles and one 
indicator based on an element-independent formulation. The element-independence of the 
new indicator is an important development because standard finite elements can now be 
used to predict error at the element-level. Use of standard finite elements is an important 
attribute for implementing this error indicator in production-level codes. 

Since the error indicator is element-level, we can compute mesh discretization error, 
element-by-element, on parallel processing machines without imposed limitations from 
multiple domain decompositions. This is an important factor in solution of large-scale 
finite element problems. 

We have defined two goals by which to conduct this research. To date, Goal 1 has been 
met and indications are that Goal 2 can likely be met. 

Comparison of the DM error indicator to analytical solutions and to well-known error 
estimators reveals that this error indicator has potential for successfully returning true 
error estimates and driving adaptive mesh refinements. 


33 



7. Topics for Future Research 

Two element-level, element-independent error indicators have been developed for use in 
various classes of problems. In addition, two adaptive mesh refinement techniques have 
been developed and used successfully. The logical completion of the current research is to 
improve both the DM error estimation and adaptive mesh refinement techniques. Com- 
pletion of the error estimation and adaptive mesh refinement research requires realization 
of the following goals: 

1. Quantitative evaluation of the current DM error indicator. 

2. Parallel adaptive mesh refinement. 

3. Application of the DM error indicator and parallel analysis procedures to more realistic 
aircraft shell structures. 

7.1. Quantitative Evaluation of the DM Error Indicator 

Numerical experiments have shown that the DM error indicator returns reasonable results 
for assessing the location of error in finite element models due to dominating modes. The 
next step will be to assess the quality of the error solution in terms of quantifying actual 
errors. For membrane-dominated problems the DM error indicator tends to underestimate 
actual error. The DM indicator has yet to be tested extensively for problems dominated 
by shear and bending. 

In order to test the DM indicator’s potential for realistic error estimation, a group of 
test problems featuring membrane bending, membrane shear, plate bending, plate shear 
and twisting will be analyzed. An initial assessment of the contribution to element strain 
energy error by each of (or combinations of) these modes can be made by considering a 
breakdown of shell element strain energies into three basic constitutive parts: membrane, 
bending and transverse shear: 

U = U m +U b + U s (7.1.1) 

where the subscripts m, b and s denote the membrane, bending and shear components, 
respectively. Computationally, this is readily accomplished for flat triangles and unwarped 
four node element by reformulating the element stiffness as 

K e = + K i + K e , (7.1.2) 

and computing the strain energy components as 1 

U = i((v e ) r K^v e + (v e ) r K£v* + (v e ) T K'v e ). (7.1.3) 

The various stiffness matrices, and hence the various strain energy components, may be 
obtained by adjusting the element constitutive model as desired for each component and 
recomputing the strain energy. Hence, by observing each of the strain energy components 
U mj U b and U s , we can acquire general knowledge about the behavior of the element and 
the source of strain energy error. It is likely, for example, that an improved quantitative 
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error estimate can be obtained by setting the shear error energy to zero in the case of thin 
shells. The decomposition (7.1.3) applies to planar shell elements in which the various 
components are uncoupled. More rigorous analysis would be necessary for warped shell 
elements in which the components are coupled. Following this basic strain energy analysis 
procedure, detailed mode-by-mode analysis will be performed. 

7.2. Parallel Implementation of Adaptive Mesh Refinement 

Parallel implementation of adaptive mesh refinement can be used to augment parallel 
solution techniques. Our interest in parallel error estimation is now extended to the parallel 
refinement of meshes. Typically, in parallel analysis, a structural model is decomposed 
into a number of substructures. This technique is a special case of the so-called domain 
decomposition (DD) method. Fig. xx shows a structure decomposed into substructures 
numbered I- VII. 

In parallel computation, each of the substructures would be assembled and solved on an 
individual processor. (The total solution might then be solved iteratively upon coupling the 
substructure solutions.) Element-level errors for the elements within a given substructure 
would also be computed on the respective substructure processor. As previously noted, the 
fact that all of the indicators developed here are element-level makes such a computation 
easily parallelizable. 

To augment the parallel solution technique with mesh adaptation, we will see that r- 
and h-adaptation schemes have severe shortcomings and that a rezoning technique would 
be more appropriate. Prior to discussion of the new technique, we briefly examine the 
problems associated with those currently implemented. 



Fig. xx Shell structure model decomposed into seven subdomain substructures. 
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7.2.1. Limitations of the r- and h-adaptation techniques 

We have shown the r- and h-adaptive refinement techniques to be effective methods for 
refining complete meshes and meshes composed of simple elements. However, both the r- 
and h-adaptive methods possess some shortcomings for mesh refinement. 

The r-adaptive mesh refinement technique requires relocation of nodes throughout the 
mesh and along external and juncture boundaries. Nodes are effectively dragged to new 
locations by error-proportionate forces from surrounding elements. In domain decompo- 
sition, juncture boundaries are often curves along which a mesh is partitioned - recall 
the example showing possible decomposition of a structural model in Fig. xx. To retain 
mesh definition and displacement compatibility (barring the use of Lagrange multipliers 
or interface elements), the juncture between sections II and IV would be a curve along 
which nodes would travel if r-adaptation were being performed. Obviously, node reloca- 
tion along the juncture would require repeated communication across the boundary which 
would violate the decomposition and incur significant processor communication overhead. 

The /i-adaptive mesh refinement technique is primarily limited to triangular elements. 
While some finite element codes employ this technique for other types of elements, 41 it 
is only done at the expense of creating elements with poor aspect or by sophisticated 
interelement constraints. Element density via h-adaptation is also limited by the size of 
the elements being split. In other words, even in areas of high error gradients, progressive 
refinements will only occur at a algebraic rate of j where h is indicative of the element size. 
Another difficulty with parallel h-adaptation is that the extensive data structure changes 
may impact load balancing unless the substructure are reconstructed periodically, which 
can be very expensive. While h-adaptation is not in direct conflict with DD analysis, its 
limitations are significant enough to consider another technique. 

7.2.2. The z-Adaptation Technique 

Based on the limitations shown for r- and h-adaptive refinement techniques, it is proposed 
to replace them with a substructure remeshing scheme known for the present as ^-adaptive 
mesh refinement. We define 2 -adaptation as “the total remeshing of a structure or sub- 
structure based on a given error distribution.” A survey of common 2 -refinement schemes 
is given by Frykestig. 42 Of the many schemes described, advancing front methods appear to 
be preferable for refinement of general surfaces. Frykestig has shown advancing front meth- 
ods work well with non-uniform rational B-spline surfaces (NURBS) 31 . These methods do 
not share many of the problems associated with other common 2 -refinement schemes (e.g., 
triangle degeneracies due to numerical round-off). A parallel implementation of rezoning 
may be performed by starting wavefront generation from each substructure boundary and 
progressing “inwards” in each. 

7.3. Analysis of Aerospace Shell Structures 

Following the evaluation and subsequent improvements to the DM error indicator, anal- 
yses of realistic aircraft-related shell structures will be performed. It is proposed that 
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general surface three dimensional test problems (perhaps derived from the CPJ bench- 
mark problem) be analyzed and compared to analytical solutions or refined solution from 
an acceptable source. 
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